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ABSTRACT 

We have analyzed new and previously published radial velocity observations of MARVELS- 1, known 
to have an ostensibly substellar companion in a ~ 6-day orbit. We find significant (~ 100 m s _1 ) 
residuals to the best-fit model for the companion, and these residuals are naively consistent with an 
interior giant planet with a P — 1.965d in a nearly perfect 3:1 period commensuribility (\Pb/P c — 3| < 
10~ 4 ). We have performed several tests for the reality of such a companion, including a dynamical 
analysis, a search for photometric variability, and a hunt for contaminating stellar spectra. We find 
many reasons to be critical of a planetary interpretation, including the fact that most of the three- 
body dynamical solutions are unstable. We find no evidence for transits, and no evidence of stellar 
photometric variability. We have discovered two apparent companions to MARVELS- 1 with adaptive 
optics imaging at Keck; both are M dwarfs, one is likely bound, and the other is likely a foreground 
object. We explore false-alarm scenarios inspired by various curiosities in the data. Ultimately, a 
line profile and bisector analysis lead us to conclude that the ~ 100 m s -1 residuals are an artifact 
of spectral contamination from a stellar companion contributing ~ 15-30% of the optical light in 
the system. We conclude that origin of this contamination is the previously detected radial velocity 
companion to MARVELS- 1, which is not, as previously reported, a brown dwarf, but in fact a G dwarf 
in a face-on orbit. 

Subject headings: stars: binaries: general — stars: individual (TYC 1240-945-1) — stars: low-mass, 
brown dwarfs 
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1. INTRODUCTION 

The measurement of precise radial velocities of stars 
has become a crucial and standard part of the exoplan- 
etary astronomer's toolkit, and its application has now 
extended far beyond its early and most successful appli- 
cation to well-studied, single, bright, chromospherically 
quiet stars. A push to apply radial velocity work to larger 
and thus fainter samples will result in the discovery rare, 
touchstone systems that will inform planet formation. 

The application to fainter stars in more crowded fields, 
as required, for instance, to follow up candidate tran- 
siting planets discovered photometrically, has necessi- 
tated attention to the effects of rare and unlikely blend 
scenarios. Since precise radial velocities can both suf- 
fer from and reveal the nature of these blended sys- 
tems, they are complementary to photometric measure- 
ments and high spatial resolution imaging using adap- 
tive optics, and thus serve a critical role in detecting, 
validating, or ruling out planet candidates in a vari- 
ety of s ituations fe.g. [Konacki et al.l [20031: iTorres et al. 
2004bt lOTJonovan et all 120071: iCollier Cameron et al. 
2007U Sant erne et al I 120121 : IBorucki et al.l 12013^ 7" Here, 
we present a case study of such an interplay in the case 
of an insidious signal seen in the MARVELS survey for 
e xoplanets. 

iLee et ail (|2011l ) announced the detection of a short- 
period "brown dwarf desert candidate" with minimum 
mass 28 ± 1.5 Mj up and period P = 5.89 d orbiting the 

F star TYC 1240-945-1 (MARVELS- 1) as the first sub- 
stellar companion discovered with the S loan Digital Sky 
Survey-in MARVELS Planet Sea r ch (lEisenstein et ahl 
[20111: IGunn et al.l 12001 IGe et al.l [200817 That work 
also reported analysis of precise radial velocities made 
with the Hobby-Eberly Telescope (HET) High Resolu- 
tion Spectrograph using an iodine cell that had residuals 
to a one-companion fit that were with surprisingly high 
(~ 100 m s _1 ). The magnitude of these residuals was 
attributed to systematic errors resulting from the pre- 
liminary nature of the Doppler pipeline, which had not 
been optimized for work at HET. 

Herein, we report followup observations with HET us- 
ing a different pipeline that confirm the high residuals 
to a single-companion fit and reveal that they are ap- 
parently consistent with an inner companion, which we 
designate MARVELS-1 c* with minimum mass of ^0.8 
Mj U p (where the asterisk indicates the provisional na- 
ture of the c component). MARVELS-1 c*, if real, would 
be unique among the known exoplanets in the size and 
proximity of its larger companion. While the period of 
this putative inner companion is ambiguous because of 
aliasing issues, the most likely solution is consistent with 
a perfect 3:1 period commensuribility (P — 1.965 d), 
indicative of a mean-motion resonance. Such a system 
would be superlative in many ways, and present a unique 
puzzle and opportunity for planet formation and system 
evolution theorists. 

We conduct an extensive examination of possible false 
positives of increasing unlikelihood, and explore the rea- 
sons why the MARVELS survey might be especially sus- 
ceptible to them. Following Sagan's maxim that extraor- 
dinary claims require extraordinary evidence, we adopt 
an attitude of healthy skepticism with respect to the 
ostensible two-companion solution. Ultimately, we find 



that the known ~6-day companion to MARVELS-1 is in, 
reality, another star in a nearly face-on orbit. 

1.1. Plan 

In Section [2] we describe our spectroscopic and imag- 
ing data of the star MARVELS-1, including its imaged 
companions and the basic stellar properties of all stars 
in the MARVELS-1 system. 

In Section [3] we describe our precise radial velocity 
measurements of MARVELS-1 with data from HET and 
Keck Observatories, our search for periodicities, and our 
best double-Keplerian orbital solution. 

In Section[4]wc describe our photometry of the system, 
which shows no variability at the 0.2 mmag at the periods 
of interest. 

In Section [5] we describe our MCMC and dynamical 
analysis of the system, which reveals only a few stable 
2-companion solutions. 

In Section[6]we discuss alternative explanations for the 
signals we see. Sections 16 . 1H6 . 51 describe a variety of false 
alarm scenarios that we can rule out or seem too unlikely 
for further consideration. Section r6.6l describes a spectral 
contamination model that describes the data even bet- 
ter than the double-Keplerian model in Section [3] Sec- 
tions [6T7] and [6]8] summarize our false alarm analysis and 
discusses the discovery of rare systems in large planet 
searches such as MARVELS. 

In Section [7] we discuss how we chose the names for the 
objects in the MARVELS-1 system, and we summarize 
our findings in our concluding Section [5J 

2. MARVELS-1 
2.1. Basic Stellar Data 

We have performed a reanalysis of the two FERO S 
spectra of MARVELS-1 presented bv ILee et ahl (|2011h . 
each of which have signal-to-noise ratios (SNR) of 340 
and resolution R ~ 48, 000, and also analyzed 2 Apache 
Point Observatory (APO) spectra with R ~ 31,500 
and each with SNR of 210. We have extracted basic 
stellar parameters using three methods: Sp ectroscopy 
Made Easy (SUE) (TValenti fc Piskunov| [l99l. StePar 
(jTabernero et al.ll2012[ and references there in) , and BPG 
(see Section 3.1.2 of IWisniewski et aT1l2012fl. T his is sim- 
ilar to the method of IWisniewski et aL (2012), but now 
includes SME as implemented at Vanderbilt as a third 
pipeline. We obtain excellent consistency in our esti- 
mates of [Fe/H] and T c g 

Averaging the six results (spectra from two sources, 
each through three pipelines), we find T c g — 6297±28K, 
log 5 = 4.22 ± 0.09, [Fe/H] = -0.13 ± 0.04, and u mic = 
1.50 ± 0.03 km s . These results h ave a significant ly 
higher gravity (a difference of 2-cr from lLee et al.ll2~011l )L 
indicating that MARVELS-1 is in fact an F9 dwarf star, 
not a sub giant as previously reported. Applying the for- 
mulae of lTorres et al.l ((2010) to these new values, we find 
M* = 1.25 ± O.O6M and R* = 1.48^;^% . 

As we will see in later sections, there is likely ~ 15-30% 
contamination in these spectra from a cooler companion, 
likely a G dwarf, that may result in additional sources of 
systematic error in these parameters. 

We have explored this possibility by analyzing the 
quality of the spectral synthesis fits generated with Spec- 
troscopy Made Easy (SME) to extract parameters of the 
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FEROS spectra. The \ 2 surfaces near minimum for these 
spectra have a parabolic character in most dimensions, 
indicating a well-behaved fit and a lack of multiple min- 
ima, as might be expected in contamination were a seri- 
ous issue. 

We quantitatively explored the effects of spectral con- 
tamination on our ability to extract precise stellar pa- 
rameter by creating a synthetic test using high-SNR 
APO spectra of HD 172051 (T cff = 5596 ± 19K,logg = 
4.56 ± 0.24) and HD 22484 (T eff = 6063 ± 19K,logg = 
4.29 ±0.16). We shifted the spectra to have the same ve- 
locity shift and scaled them so that the former (cooler) 
star contaminated the latter with a contribution of 30%. 
(We chose this level of contamination because it approx- 
imates the flux ratio expected if these particular stars 
were at the same distance, and is at the upper end of the 
contamination we expect in MARVELS-1). 

We then used the StePar code to estimate the stellar 
parameters of the blended spectrum and found T e s = 
5987 ± 23K, log# = 4.44 ± 0.19. These values are, not 
surprisingly, intermediate to the correct values of the con- 
tributing spectra. In the case of T e ff the value from the 
blended spectrum is closer to the brighter star, and in the 
case of log g the value from the blended spectrum is con- 
sistent with that of the brighter star within the quoted 
uncertainties. 

This test is not conclusive because the "secondary" star 
is more metal poor than the primary by 0.2 dex ([Fe/H] = 
-0.29±0.02 vs. -0.07±-0.02), which may effect the final 
result (which yielded [Fe/H] = -0.11 ±0.05), but it does 
suggest that the values we derive here for MARVELS-1 
are not badly affected by a contaminating spectrum. Ad- 
equate caution is therefore needed when using these pa- 
rameters, which were calculated assuming an unblended 
spectrum. 

2.2. Companions Revealed by AO Imaging 
2.2.1. Data acquisition and reduction 

We acquired images of MARVELS-1 using NIRC2 (in- 
strument PI: Keith Matthews) with the Keck II adaptive 
optics (AO) system ()Wizinowich et al.|[2000h on UT 24 
June 2011. Our initial set of observations consisted of 
dithered images taken with the K' (A c = 2.12/im) filter. 
Using the narrow camera setting, to provide fine spa- 
tial sampling of the NIRC2 point-spread function (PSF), 
we obtained nine frames with 4s of on-source integra- 
tion time each. A preliminary inspection of raw data 
revealed the presence of an additional point source lo- 
cated 0'.'900 east of the primary. Upon noticing the can- 
didate companion, we obtained images with the J-filtcr 
(A c = 1.25yum) to provide complementary photometry, 
and detected the companion at 0'.'897 separation, consis- 
tent with the K' measurement. 

We processed the images using standard techniques to 
replace hot-pixels, flat-field the array, and subtract ther- 
mal background radiation. After aligning and coadding 
the frames, we noticed a third object in the field, to the 
north of the primary star and separated by only 0'.'152 
in K' band (and ~0'.'18 in J band, though we adopt the 
K 1 band value because the PSF subtraction is cleaner in 
those images). 

Figure Q] shows the image of MARVELS-1 in the K' 
filter. Arrows indicate the location of the 0'.'9 and 0'.'15 



companions. Inspecting the multi-color data set, we find 
that neither object's position moves as a function of 
wavelength, demonstrating that they are not speckles. 
Further, the PSF of each source is qualitatively similar 
to that of the primary star. 

In the rest of this manuscript, we will refer to 
the primary star as MARVELS-1 A, the 0'.'9 compo- 
nent as MARVELS-1 C, and the 0'.'15 component as 
MARVELS-1 B (see Section [7] for an explanation). 

We performed aperture photometry to further charac- 
terize each point source using the publicly available tool 
Starfinder (jDiolaiti et al.ll2000h . Often used for studying 
the Galactic center and other crowded fields, Starfinder 
is optimized for AO photometry for which numerous 
sources are spatially blended. By self-consistently fitting 
each PSF, using an empirical iterative algorithm that 
minimizes residuals, we were able to account for con- 
taminating light from the much brighter primary star to 
measure an accurate relative brightness, angular separa- 
tion, and position angle of each candidate in each filter. 
Our astrometric and photometric measurements for both 
companions relative to the primary are summarized in 
Table [Q MARVELS-1 C has colors comparable to an 
M9V dwarf (J - K = 1.245±0.038), and MARVELS-1 B 
has colors (J — K = 0.804±0.071) consistent wi th spec- 
tral type K7V-M2V (H raus fc Hillenbrandl[2007l) . 

We have in terpolated the absolu t e mag nitude values 
in Table 5 of iKraus fc Hillenbrandl (|2007n for the spec- 
tral type of MARVELS-1 A (F9V, Section EH, and 
derive Mr- =3.11 and Mj =3.42. From the apparent 
magnitudes of MARVELS-1 (K = 9.032 ± 0.017 and 
J = 9.395±0.018 according to lLee et al.ll2011f ). we derive 
a rough distance of d = 150 ± 30 pc from spectroscopic 
parallax. 

MARVELS-1 C is too bright to be an M9V dwarf at 
the same distance as MARVELS-1 A, it is therefore likely 
an unbound foreground star. The colors and absolute 
magnitude of MARVELS-1 B are consistent with the as- 
sumption that MARVELS-1 A and B are at the same 
distance. 

We obtain a consistent mass esti mate for MARVELS-1 
B of O.62M from the relations of iDelfosse et al.1 ()2000[ ) 
in both J and X-band (using Mj — 5.95 and Mr- = 
5.20). 

MARVELS-1 has a proper motion of alpha = 33.6, 8 = 
—5.8 mas yr _1 . With more than a one year time baseline 
from our first epoch observations, future precision astro- 
metric measurements using NI RC2, which has a plate 
scale of 9.963 ± 0.006 mas (jGhez et al.ll2008l ) will be suf- 
ficient to determine whether the two visual companions 
share a physical association with the primary. The angu- 
lar proximity and the consistency of the MARVELS-1 B 
magnitudes and colors with being a dwarf at the distance 
of MARVELS-1 A argues that they are likely bound. Ir- 
respective of whether these sources represent a chance 
alignment along the line of sight, their flux contribution 
must be accounted for when analyzing precise Doppler 
data. 

3. RV OBSERVATIONS 

3.1. HET Confirmation Data 

iLee et all (|201 If) reported nine precise radial veloc- 
ity measurements made with the Hobby-Eberly Tele- 
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Our new observations were made in late 2010 and early 

TABLE 1 2011, sufficiently long after the original observations to 
Relative Astrometry and Photometry for the AO-resolved . , ° r , , , , ., , 
components of MARVELS- 1 with respect to the A component provide some coverage of previously unexplored orbital 
phases. This system is particularly challenging to ob- 



component 



AK' 
mag. 



AJ 
mag. 



Separation 
arcsec. 



PA 



NIRC2 , K 

June 24, P 2011 UT 



serve at HET because the companion periods are so close 
to integer numbers of sidereal days, and the fixed eleva- 
2.088±0.035 2.529±0.056 0.153±0.0017 2.6±0.33 tlon of HET requires that all observations occur within 
C 3.429±0.02l 4.311±0.018 0.8998±0.0003 86.25±0.02~ 1 hour of two sidereal times (corresponding to the 

rising and setting tracks of MARVELS-1). This con- 
straint made acquiring good phase coverage problem- 
atic, and complicated our ability to rule out competing 
and qualitatively different orbital solutions, and to deter- 
mine orbital eccentricities. Nonetheless, we were able to 
adaptively schedule o ur observations in the HET queue 
(|Shetrone et al.ll2~007l) so as to optimally explore the 2 d 
period and eliminate all but one among many competing 
orbital solutions near 2 d. 

We also obtained a template and 14 velocity mea- 
surements from Keck using the slit-fed HIRES spectro- 
graph in August and September 2011. These observa- 
tions were reduced us ing the usual methods of the Cali- 
fornia Planet Survey (IHoward et all [2010; Joh nson et al.l 
I2010t I Wright et al.l I201ll and references therein) . Due 
to the higher signal to noise ratio of these observations 
compared to the HET data, the corresponding internal 
errors on the radial velocity measurements are signifi- 
cantly smaller. 



+ 
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Fig. 1. — K' band adaptive optics image from Keck/NIRC2 show- 
ing MARVELS-1 B and MARVELS-1 C. 

scope (iRamsev et al.l I1998D High Resolution Spectro- 
graph (HRS; iTulJ U9980 in late 2009 with Director's 
Discretionary Time as part of the radial velocity con- 
firmation program to the discovery of MARVELS- lb 
made with the MARVELS instrument. These nine ob- 
servations were all acquired in 2 009 December with an 
iodine absorption cell calibrant (|Butler et al.l I1996D at 
i?- 60, 000. 

These nine observations were reduced with a prelimi- 
nary version of a precise Doppler pipeline (provided by 
Debra Fischer). This pipeline included an instrumental 
profile model more appropriate for the slit-fed Hamilton 
spectrograph at Lick Observatory and other components 
but was not optimized for work on HRS. 

These nine points had significantly higher native pre- 
cision tha n the SMARTS and MARVELS velocities re- 
ported by iLee et al.l (| 2 1 If ) . However these velocity er- 
rors had to be scaled by a factor of ~ 15 to match the 
observed ~ 100 m s _1 scatter about the single-object 
Keplerian fit to the data (though with a 6-component 
Keplerian fit to only 9 measurements, the residual scat - 
ter was difficult to quantify precisely). ILee et all (|2011D 
attributed these high residuals to systematic errors in 
their less-than-optimized Doppler pipeline. 

3.2. New Data 

The large residuals in the HET confirmation result did 
not appear completely random and exhibited significant 
power at 0.5 d _1 . This and a desire to diagnose the 
remaining systematic errors in the Doppler pipeline led 
us to request Director's Discretionary Time on HET to 
acquire 21 additional epochs on MARVELS-1 at higher 
signal-to-noise ratios. These new observations were made 
with the same spectrograph settings as the original nine 
observations. 



3.3. Data Analysis: Template Analysis, Raw 
Reduction, and Doppler Pipelines 

We performed raw reduction of th e echellograms from 
HRS u sing the REDUCE package of IPiskunov fc Valentil 
( 2002) . This package required some modifications to the 
slit function description to accommodate the flat-topped 
n ature of fiber-fed sp ectra. 

iWang et al.l (|2012l) describes our Doppler pipeline. In 
brief, the code base developed by John Asher John- 
son and the California Planet Survey (e.g. IHoward et all 
l200l I2011bl lal) was modified for the HRS at HET This 
code is based on the principles of lButler et al.l (|1996l ). We 
have demonstrated its precision on HRS spectra with the 
RV standard star (a Dra), for which we achieve an r.m.s. 
scatt er below 3ms" 1 (|Wang k. Wrightl2011l:IWang et al.l 
I2012T) . 

Our HET observations of MARVELS-1 arc not strictly 
analogous to our a Dra observations, however, in that 
we have a lower signal-to-noise ratio in both the template 
and radial velocity observations. To put a very conserva- 
tive upper limit on the systematic errors inherent in our 
Doppler pipeline on fainter targets, we analyzed three 
other MARVELS targets, TYC 1194-144-1, TYC 3413- 
2471-1, and TYC 3410-1406-1, in a manner identical to 
MARVELS-1. 

We stress that these targets were observed at signifi- 
cantly lower signal-to-noise ratios than MARVELS-1 as 
exploratory science and so their RVs should have consid- 
erably higher random and systematic noise. These three 
stars represents our entire HET sample of non-binary 
MARVELS targets, at present; we have no reason to con- 
sider these to be RV stable stars, since these are, to our 
knowledge, the first precise radial velocity measurements 
ever obtained on these targets. 

We present these three radial velocity targets in Fig. [2j 
The noisiest velocities here are for TYC 1194-144-1, with 
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an r.m.s. dispersion of 25 m s , the best are for TYC 
3413-2471-1, at 9.5 m s" 1 . Our MARVELS- 1 observa- 
tions were made at higher signal-to-noise ratios and we 
have no reason to suspect MARVELS-1 has levels of jitter 
this large; we thus expect to achieve significantly better 
precision on MARVELS-1, closer to our performance on 
a Dra. 

We obtained an iodine-free stellar template at 
R=120,000 for this target, but the signal-to-noise ratio 
was low due to the faintness of the star and the slit losses 
at this high spectral resolution at HRS. We opted there- 
fore to use the deconvolved stellar template obtained 
from Keck observatory. This template produced supe- 
rior interal errors in the HET velocities. 

All of our HET and Keck velocities appear in Tables [2] 
& [3] The reported uncertainties are "internal" errors 
measured by the consistency with which different por- 
tions of the spe ctrum ( "chunks" ) report the radial veloc- 
ity of the star (|Butfer et al.lll996l ). 

3.4. Keplerian fit 

We first performed a standard double-Keplerian fit, 
with no accounting for dynamical (Newtonian) interac- 
tions or other sources of non-Keplerian signal. 

The b component of the MARVELS-1 system has such 
a large amplitude that its orbital parameters are insen- 
sitive to the fit for the c* com ponent, and rema in largely 
unchanged from the values in ILee et all (|2011l ) . Fits for 
this system are complicated by the fixed zenith angle of 
the HET (see § I3.2p which makes it difficult to acquire 
good phase coverage for this object, and produce signifi- 
cant side lobes in the periodograms of residuals to a one 
companion fit, as shown in Fig. [3] 

The periodogram shows two dominant peaks, near 2 d 
(0.5 d _1 ) and 0.66 d (1.5 d _1 ), which have similar am- 
plitude and each with substantial numbers of sidelobes. 
These two peaks are aliases of each other (their frequen- 
cies differ by exactly 1 sidereal day -1 ). As we describe 
below, only the 1.965 day peak corresponds to orbital 
solutions with residuals near our expectations; all of the 
other peaks (including the 0.66 d peak) have best fits 
with at least 10 m s _1 higher levels of residual scatter. 
Given that the 2 d peak has a perfect period commensu- 
ribility with the primary signal and the 0.66 d does not, 
the 2 d peak is more likely to be the "real" signal. 

Regardless of its true period, the new observations ap- 
pear to strongly confirm the presence of a second peri- 
odic signal of apparent semiamplitude ~ 100 m s _1 for 
MARVELS-1. W e have fit all velocity poi nts using the 
RVLIN package of I Wright fc Howard! (|2009l ). checking for 
a suite of orbital periods near 2 d and 0.66 d to explore 
all of the peaks shown in Fig. [3] We assumed 9 m s _1 
of jitter, chosen because it produced a xl near 1 (the fit 
details are not strongly sensitive to the amount of jitter 
assumed, and we performed a more robust Markov-chain 
Monte Carlo (MCMC) calculation with jitter as a param- 
eter, as well, described in § 15. ip . The parameters of this 
best 2-companion fit are shown in TableSJ We calculated 
parameter uncertainties for this fit u sing bootstrappin g 
methods with the BOOTTRAN routine ()Wang et alj2 012). 

Figure H] shows the measured velocities after subtract- 
ing the solution for the b component from the best 2- 
companion fit. The temporal baseline of the velocities 



TABLE 2 

HET Radial Velocities for MARVELS-1 
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64 


14 


6 


15483.7415 


2159 


31 


13 


5 


15483.7528 


2142 


75 


11 


7 


15484.9464 


-489 


47 


13 


5 


15484.9577 


-507 


08 


12 


7 


15485.7397 


-2353 


73 


Id 


4 


15485.7510 


-2359 


49 


17 


4 


15497.7059 


-2631 


09 


18 


7 


15498.7192 


-2078 


22 


15 


5 


15498.9215 


-1623 


45 


20 


6 


15499.7351 


359 


36 


16 


1 


15500.6947 


2311 


16 


18 


1 


15500.9073 


2468 


96 


17 


6 


15501.6961 


1699 


77 


18 


1 


15507.6863 


1526 


22 


14 


1 


15510.6632 


-1762 


64 


15 





15510.6715 


-1711 


09 


15 


5 


15522.6391 


-1285 


58 


18 


2 


15522.6469 


-1283 


72 


16 


5 


15522.6547 


-1258 


66 


17 


9 


15522.8492 


-777 


53 


17 


7 


15522.8570 


-728 


73 


15 


8 


15523.6459 


1120 


41 


13 





15524.6390 


2472 


26 


17 


2 


15524.8380 


2341 


07 


18 


8 


15527.6284 


-2843 


01 


15 


2 


15531.6080 


708 


58 


17 


2 


15576.7008 


1109 


81 


21 


7 


15577.7000 


2483 


02 


16 


9 



span over a full year, giving coverage of most phases of 
the orbit at this best fit period. The r.m.s. to the fit (the 
standard deviation of the residuals, with no adjustment 
for the number of model parameters) is 14.7 m s . 

Under the assumption that the c* component is real, 
we would attribute this 9 m s -1 of jitter needed to 
achieve xt near 1 to strong gravitational perturbations 
fsee SIOl). 

4. PHOTOMETRY 

The transit sea rch photometry of Lee et al.l ([201 If) us- 
ing KELT data (jPepper et all 120071: iSiverd et all 120091) 
rules out photometric variability on the period of either 
RV signal at the sub-mmag level; their Figure 3 illus- 
trates that the star is photometrically constant. 

We have reanalyzed these data, and find that there is 
no evidence for periodic variability, and the limits for 
P=l to 10 days are quite stringent. The strongest peak 
in the periodogram in that interval has P = 1.0218d 
with amplitude 0.87 ± 1.8 mmag, and we associate this 
peak with diurnal effects. At the periods of the signals 
of MARVELS-1 b and c*, the best-fit photometric am- 
plitude is < 0.2 mmag. We show the KELT data phased 
a t these per i ods in Figure [Sj 

ILee et ail ()201lD used these data to search for transits 
of the b component. They were able to rule out transits 
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Fig. 2. — Radial velocities for three other MARVELS targets of similar apparent visual magnitude to MARVELS- 1, but observed at 
lower signal-to-noise ratios and analyzed using noisier template observations. These observations represent a conservative upper limit to 
any systematic errors in our Doppler pipeline as applied in practice to typical MARVELS targets. We have no a priori reason to believe 
that these stars are inherently RV stable, and so the r.m.s. scatter of these radial velocities (listed in each pane in m s -1 ) represent upper 
limits to our precision on these targets. We expect to achieve significantly better precision on MARVELS-1, closer to our performance on 
a Dra, for which pipeline achieves 3 m s _1 . 




Period Period 

Fig. 3. — Left: Periodogram of the residuals to a one-companion fit to the HET RV data. There is substantial power near 2 days, but 
severe aliasing due to the near-integer period of the c* component produces a second peak near 0.66 days and large sidclobes to the tallest 
peaks. The 2-companion Keplerian fit with the smallest residuals has P c = 1.965 d. Right: Detail near 1.965 d. 



with a depth of of ~ 0.2% or greater, corresponding to 
companions of Jupiter radii or larger for the original es- 
timate of the radius of the primary. If the c* component 
exists, its a priori transit probability is ~ 20%, and for 
our new stellar radius, a Jupiter radius companion at the 
period of the c* component would produce a transit with 
a duration of ~ 3 hours and a depth of ~ 0.5%. For a 
strictly periodic transit, we can robustly exclude such a 



transit signal. However, we expect any transits of c* to 
exhibit large TTVs with amplitudes of up to 10 hr due to 
dynamical interactions with b (see Section 15. 3p . There- 
fore, when phased to a single period, the transit will get 
smeared out, thereby lowering the signal-to- noise ratio 
with which we detect or exclude the signal. 

Nevertheless, we expect that such transits would prob- 
ably still be detectable in the binned curves, even given 
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Fig. 4. — Solid black disks are HET velocities, crosses are Keck velocities, and small gray symbols are repeated points. Left: Radial 
velocities phased at the period of the b component. On this scale, the effects of the c* component are just visible. The internal (random) 
errors are smaller than the large plotted points. Right: Residuals phased to the period of MARVELS- 1 b given in Table [4] showing three 
consecutive periods (only the gray points are repeated). Deviations from a Keplerian curve are large, and appear to be a consistent function 
of phase of the b component. Note also that the zero crossings occur at the simultaneously with the zero crossings of the b component. We 
phase both diagrams to have phase zero at the velocity zero with positive slope of the primary component. 




Fig. 5. — KELT photometry of MARVELS-1 phased at the periods of the b (left) and c* (right) components. Gray points show the 
individual measurements and the black points show the data binned. The binned data show that the star is photometrically stable at the 
sub-mmag level, with no evidence of either transits or stellar pulsations. 



TABLE 3 

Keck Radial Velocities for MARVELS-1 



BJD 



time 
-2440000 (UTC) 



Velocity 
m s — 1 



Uncertainty 
m s _1 



15769. 
15782. 
15783. 
15789. 
15791. 
15793. 
15796. 
15797. 
15798. 
15808. 
15809. 
15810. 
15811. 
15812. 



0930 
1024 
1275 
1125 
1307 
1094 
1305 
1316 
0922 
0954 
0594 
0965 
0929 
1337 



-2349.90 
-680.13 
1666.41 
1886.25 
838.79 
-2244.26 
2644.81 
599.29 
-1702.91 
2335.50 
266.61 
-2077.35 
-1821.60 
643.32 



5.6 
4.3 
5.2 
5.5 
3.3 
5.4 
5.9 
2.8 
5.5 
6.5 
3.1 
5.6 
4.6 
3.6 



the large possible transit timing variations: we can ap- 
proximate the signal as being diluted by a factor of the 
duration/TTV amplitude, which is ~ ( Sh/lOh) * 5% ~ 
0.15%, a level similar to that for which ILee et al.l ((201 lh 
ruled out transits of the b component. 

5. MCMC AND DYNAMICAL ANALYSIS 

The novelty of this system demands a thorough anal- 
ysis of the orbital parameters of the components and 
their dynamics. To this end, we explored the space of 
two-companion orbital fits to the RV data with a combi- 
nation of MCMC analysis of the a posteriori Keplerian 
orbital elements and with a dynamical analysis of some 
of the families of solutions. This MCMC analysis was 
performed at a relatively early stage of the investigation, 
before the Keck data was available. As such it is based 
on HET data alone. The MCMC analysis was not re- 
newed following the Keck data observations, since those 
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TABLE 4 

Formal best-fit Keplerian Orbital Elements 
for substellar companions in the 
MARVELS-1 System. Note that the 
parameters of our best model for. the system 
appear in Table 7, not here. 



Property 


b 


c* 


Per (d) 


5.895394(63) 


1.96492(29) 


T (JD-2440000) 


15498.059(91) 15498.770(91) 


e 


0.0160(16) 


0.134(56) 


w (°) 


179.7(5.5) 


306(19) 


K (m s- 1 ) 


2572.6(4.0) 


104.8(4.4) 


msini (M Jup ) 


27.6(1.5) 


0.764(42) 


a (AU) 


0.0702(19) 


0.03351(91) 


r.m.s. (m s" 1 ) 


14.7 






1.01 




jitter (m s -1 ) 


9 




A'obs 


30 




Note. — This fit 


is dynamically unstable. For 



succinctness, we express uncertainties using paren- 
thetical notation, where the least significant digit 
of the uncertainty, in parentheses, and that of the 
quantity are to be understood to have the same 
place value. Thus, "0.100(20)" indicates "0.100 ± 
0.020", "1.0(2.0)" indicates "1.0±2.0", and "1(20)" 
indicates "1 ± 20" . 

data revealed that the c* signal was likely spurious. 

The period ambiguity from the aliasing of the pe- 
riod with the telescope observing constraints prevent our 
MCMC code from settling on a unique period for the c* 
component, consistent with the heights of the sidelobes 
in Fig.[3J We choose to focus on the 1.965 d period, which 
has the lowest residuals and is favored by the data. This 
inner period is consistent with a perfect 3:1 PC (the best 
fit Keplerian solution has \Pb/P c — 3| ~ 10~ 4 ). 

In summary, we find that most of the solutions in the 
MCMC chains are unstable on short timescales, but there 
do exist stable solutions. If the solutions we explore 
here describe an actual qualitative behavior of the sys- 
tem, then this system would exhibit significant perturba- 
tions, which would imply that the osculating Keplerian 
orbital elements would vary strongly and detectably on 
timescales of years, and that planetary transits might be 
intermittent as the inclinations of the companions oscil- 
late. Were transits to occur, they might exhibit large 
(~ 10 hour) timing variations. 

5.1. MCMC Methodology 

We analyzed the radial veloci ty me a surem ents using 
a Baye sian framework following Ford (2005) and iFordl 
(|2006aD . We assume priors that are uniform in eccentric- 
ity, argument of pericenter, mean anomaly at epoch, ve- 
locity zero-point and logarithm of the orbital period. For 
the velocity amplitude (K) and jitter (cj), we adopted 
a prior of the form p(x) = (x + x ) _1 [log(l + x/xo)}^ 1 , 
with K Q — Oj, = 1 m s _1 , i.e., high values are pena l- 
ized (for a discussion of priors, see lFord k. Gregorvl 2007). 
The likelihood for radial velocity terms assumes that each 
radial velocity observation is independent and normally 
distributed about the true radial velocity with a variance 
of af + o~j, where <7j is the published measurement un- 
certainty, and a j is a jitter parameter that accounts for 
additional scatter due to stellar variability, instrumental 
errors and/or inaccuracies in the model (i.e., neglecting 



perturbations or additional, low amplitude planet sig- 
nals) . 

We used an MCMC method based upon Keplerian or- 
bit fitting to calculate a sample from the posterior dis- 
tri bution (IPordl [2006b). We use the algorithm described 
in IFordl (|2005[ ) to adaptively determine the appropriate 
step size to employ in our Markov chains, initializing 
the algorithm with system parameters drawn from the 
best fits found using the frequentist approach described 
in section 4. We calculated 8 Markov chains per MCMC 
realization, each with ~ 2 x 10 s states, and discarded 
the first half of the chain s. We calculated Gelman-Rubin 
(|Gelman &: Rub in 1992) test statistics for each model pa- 
rameter and several ancillary variables and found no indi- 
cations of non-convergence amongst the individual chains 
(but see section 15.21 for discussion of the differences be- 
tween individual fits). 

Following Keplerian fitting procedure, we attempted 
n-body DEMCMC fitting of the system (using the 
method described in Ite r Braakl2006HPavne fc Fordl l2011t 
Uohnson et al.ll201l|) . but found that either, (i) the total 
number of observations is currently too small to constrain 
the parameter space sufficiently for this numerically in- 
tensive mechanism to converge in a reasonable time, or 
(ii) no self-consist planetary fit is possible for this system. 

To ensure that at least some fraction of the Keplerian 
fits are stable in the strongly-interacting regime expected 
of objects that are so close together and so massive, we 
took the results of the MCMC fits and i njected those sys- 
tems into the Mercury n-body package (|Chamberslll999| ) 
and integrated them forward for ~ 10 5 yrs (~ 10 7 or- 
bits). This exercise allowed us to distinguish the long- 
term stable fits to the data from those systems which 
happen to fit well but become unstable on very short 
timescales. We define "unstable" in a loose manner: we 
reject any systems which undergo a collision or change 
either of the planet's semi-major axes by > 50%. This 
approach allowed significant variation in elements to oc- 
cur, but rejects systems which qualitatively change their 
architecture. 

We assumed that all systems are coplanar and edge-on 
for the sake of this analysis, hence all of the masses used 
in our n-body analyses are minimum masses. 

5.2. MCMC Results 

We performed a number of independent MCMC sim- 
ulations, all of which appeared to converge individually. 
However, the point of convergence often varied in dif- 
ferent runs, which manifested itself in different sets of 
solutions having slightly different periods for the c* com- 
ponent (i.e. the chains converged on local minima at 
the posterior probability distribution) essentially corre- 
sponding to the sidelobes for the inner period evident in 
Fig. 13 

For the rest of this analysis we concentrate on those 
chains that converged on the local minimum character- 
ized by the lowest values of the jitter (~5ms _1 ). These 
chains had an associated period for the c* component of 
P c ~ 1.965 d. The chains which converged to alternative 
local minima had significantly larger jitters (~ 10 — 50 
m s _1 ). The prior for jitter is intentionally chosen to 
have significant support even at realistically large jit- 
ters (which are highly unlikely due to the absence of 
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TABLE 5 

Orbital elements at JD 2,455,200 (noon 3 Jan 
2010 UT) for the substellar companions 
resultin g fr om the keplerian mcmc 
analysis of 35.21 in which jitter is included as 
a free parameter. note that the parameters 
of our best model for the system appear in 
Table 7, not here. 



Quantity 


Mean 


+ 1<T 


-lo- 


P b (d) 


5.89537 


0.00007 


-0.00007 


K b (m s" 1 ) 


2572.4 


4.3 


-4.4 


e b 


0.0158 


0.0017 


-0.0017 


"b (°) 


179.4 


6.0 


-6.1 


M h (°) 


202.3 


6.0 


-6.0 


m b sini b (M Jup ) 


28.01 


0.05 


-0.05 


Pc (d) 


1.9649 


0.0003 


-0.0003 


K c (m s- 1 ) 


103.1 


5.0 


-4.7 


e c 


0.101 


0.056 


-0.056 


^ (°) 


279.1 


27.9 


-40.1 


M c (°) 


113.2 


28.1 


-29.6 


m b sini b (M Jup ) 


0.778 


0.036 


-0.034 


Log Jitter 


0.185 


1.58 


-2.41 


x 2 


3.711 


0.527 


-0.52 



Note. — Many of the solutions spanned by 
these parameters are dynamically unstable on short 
timescales. 



strong stellar activity indicators in this star) so as to 
minimize the risk of overlooking interesting structure in 
the goodness-of-fit surface. 

The results of our Keplerian MCMC fitting are pre- 
sented in Fig [6] and summarized in Table [5j This 
analysis favors P c = 1-964931°;°°^° d and P b = 

5.895374t°;oooo66 d > implying a period ratio that is 
consistent with exactly 3:1, within the uncertainties 
(Pb/Pc = 3.00029721+°;^). The MCMC fits con- 
strain the outer object to have a low eccentricity e& = 
0.0161qqq2, but allows a large spread in eccentricity for 
the c* component, e c = O.lOt ,; ^. 

5.3. Dynamical Analysis 

As outlined in H5.ll the results of the Keplerian MCMC 
analysis in £15.21 were used as inp ut to the Mercury n- 
body integrator (|Chambersl I1999T) to integrate the sys- 
tems for 10 5 years. "Stepsize chaos" can be a source 
of error in numerical integrations, but it is negligible 
when the timestep At is smaller than the shortest phys- 
ical timescakjjn_the^y^tem by at least a factor of 10 or 
20 (|Rauch fc Holmari ll999'). Each integration was car- 
ried out using a fixed timestep of At = 0.01P C , where 
P c is the period of the inner planet. This exercise dra- 
matically reduces the number of suitable orbital fits, as 
> 99% of the solutions become unstable (using the ap- 
proximate definition in £I5.1|) over the course of the 10 5 
yr simulation. 

The stable systems favor lower periods and eccentric- 
ities for the c* component, as shown in the first col- 
umn of Table [6] (we stress that the small number of sys- 
tems (~ 20) that remained stable means that the val- 
ues quoted suffer from small- number statistics). While 
the best fit values for the period and eccentricity of the 
outer object are essentially unchanged, the best fit val- 
ues for the c* component become P c = 1.9644 + .° ) ; , ) )2 d 



and e c = 0.017_q'q}4, respectively. This means that the 
period ratio is increases slightly to ~ 3.001, but, more 
significantly, a much more circular orbit for the c* com- 
ponent is strongly preferred (at least at the epoch of first 
observation — see below for further discussion). 

The results of the long-term stability analyses enumer- 
ated in the first column of Table [5] should not be taken 
to be our final determination of the system parameters. 
The low fraction of the systems from the Keplerian fit 
which were stable indicates that while it is reassuring 
that there exist stable solutions within the neighborhood 
of the Keplerian fits, further work must be done to pro- 
vide a reasonable statistical sample of such orbits. The 
orbits are not Keplerian over the duration of the fits (see 
below for further discussion), so a rigorous DEMCMC 
n-body fitting needs to be undertaken to generate such 
a sam ple of stable orbits (as outlined in Uohnson et al.1 
2011). As noted in our attempts at such an anal- 
ysis were unsuccessful, as the routines did not converge. 
This could either be due to a lack of a sufficient number 
of high cadence observations, or may simply be due to 
it being impossible for a stable planetary system to give 
rise to the observed RV signal. 

Fig [JJ shows the detailed evolution over the course of 
10 5 yrs of various orbital elements in the system for one 
of the stable solutions. In addition to the period ratio, 
pericenters, semi-major axes and apocenters, we also plot 
the resonant argument 0\ , 02 and #3 (Murray & Dermott 
2000), defined af3 

9\ = 3Aj — A c — 2w c 

62 = 3A b - A c - (lj c + uj b ) (1) 
03 = 3Ah — A c — 2wf, 

One can see that the system is in a resonant configura- 
tion, with 02 (gray asterisks) librating about 0° with an 
amplitude of ~ 90°. One expects objects in resonance to 
have semi-m ajor axes which oscillate slightly as a func- 
tion of time (|Murrav fc Dermottlll999l ). as demonstrated 
in the middle panel of Fig [7] Interestingly, in the case 
of MARVELS- 1 this leads to period oscillations of ~ 3% 
over the course of ~ 200 days (< 40 orbits of the outer 
body). Similarly, the eccentricity of the c* component 
in this sample simulation changes between ~ 0.01 and 
- 0.1 in- 100 d. 

If this example simulation were typical of the real 
system, then such oscillations in period and eccentric- 
ity would be the largest and fastest seen to date in any 
exoplanetary system. This result would have several im- 
plications: 

i For systems which interact stro ngl y on short 
timescales. lLaughlin fc Chambers! (|2001l) and 
iRivera fc Lissauerl (|2001[ ) demonstrated that Keple- 
rian orbital fits are insufficient and that interaction 
terms need to be accounted for to fit the radial 
velocity data. While this is patently the case here for 
MARVELS-1, in this situation the interactions are 
so strong and orbital perturbations occur on such 
short timescales that one can only report the range, 
manner and timescale over which the quantity varies, 

19 where A is the usual mean longitude and we have taken the 
longitude of the ascending node Q = for simplicity. 
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Fig. 6. — The results from the Keplerian MCMC fitting algorithm which converged to low jitter solutions. The top row contains (Pi,, P c ) 
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the mode. We find that the period of the b component is very tightly constrained to be P;, = 5.895374 
is less tightly constrained, lying in the range P c = 1.96492 



+0.00029 



+0.000066 
-0.000067 



d, while the c* component 



in a manner similar to that undertake n for the proper 
orbit al elements of asteroids (e.g. iKnezevic et al.l 
l200l . 

ii The high masses and tight orbits of the bodies in 
TYC 1240 may provide an even more exquisite probe 
of a dynamicall y interacting system than the famed 
GJ 876 system rtMarcv fc Benit3ll989HDelfosse et al.l 
1998; Ri vera et al.ll2010D : if a sufficient number of suf- 



d, where the uncertainties indicate ~ l-cr confidence levels. 

ficiently precise RV observations can be taken with a 
sufficiently high cadence, then the dynamical inter- 
actions between the companions could reveal their 
inclinations and true masses. The advantage of the 
MARVEL S-l system is that nearly a full cycle of the 
orbital element variations takes place within a single 
observing season, significantly boosting observational 
plausibility. 
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TABLE 6 

Orbital elements for the substellar companions resulting from integrations 
of the long-term stable results arising from the keplerian mcmc analysis 
of £15.21 Note that the parameters of our best model for the system appear 

in Table 7, not here. 



Quantity 



Mean Value at Epoch 1 1-Year Mean 2 1-Year Min 2 1-Year Max 2 



Pb (d) 
a b (AU) 
K b (m s" 1 ) 


t>.»y;xS4_ 00005 

n 071 977R+0. 0000005 
U.U(1Z( (D_ 0000004 

2572.4 + || 


5.944(1) 
0.071332(8) 


5.927(1) 
0.071200(10) 


5.961(1) 
0.071468(9) 


e-b 


0.0168+^ 


0.020(2) 


0.016(2) 


0.024(2) 




145.1+g 3 


176.5(6.3) 


162.8(8.4) 


189.1(5.9) 



Pc (d) 


1 9644+ auuu2 
i.»o<±^_ 0002 


1.943(2) 


1.904(3) 


1.981(4) 


a c (AU) 




0.033846(29) 


0.033397(39) 


0.034288(48) 


K c (m s" 1 ) 


103.2111 








e c 


o m7+ 0013 

u - ui '-0.014 


0.058(7) 


0.0006(3) 


0.14(1) 




252 4+ 99 - 1 
zdz -^-192.8 


143.0(36.5) 


0.0(1) 


360.0(1) 



Note. — These elements span our favored dynamical solutions, however we cannot 
yet conclusively describe the true qualitative behavior of the system. In addition, the 
perturbations in these solutions make these Keplerian osculating elements an inade- 
quate description of the system. See Table [4] for a note on error notation. 

1 Results for JD-2,455,200, averaged over the 20 results from the MCMC fitting of 
Table [5] that were also long-term stable. 

2 Extracted from detailed 1-year simulations initialized from the 20 individual starting 
conditions which underly the figures presented here in the "Mean Value at Epoch" 
column, with data output at hourly intervals and means, minimums, and maximums 
calculated. 



iii Such strong interactions and associated large changes 
in orbital elements will likely drive significant transit 
timing variations (TTVs) in any hypothetical transit 
of the c* component. 

iv Our calculations were specifically for coplanar sys- 
tems; however, if the companions are mutually in- 
clined, then the precession of the orbit of the c* com- 
ponent could be sufficiently strong and rapid that the 
system could quickly pass through episodes of tran- 
sit and non-transit, further complicating any transit 
search. 

We address point (i) above in Table [51 We have taken 
all 20 stable systems from our long-term integrations 
and extracted the osculating elements at the epoch of 
the first observation. We report the mean and standard 
deviation of these parameters across the 20 solutions in 
the first data column. To capture how they vary with 
time, we have also calculated the mean and extremum 
values of these elements over the course of each of the 20 
simulations, and report the mean and standard deviation 
of these quantities across the 20 simulations in the other 
three data columns. 

For instance, over the period of the n-body calculations 
of the stable solutions the eccentricity of the c* compo- 
nent has a mean value of e c = 0.058 and minimum and 
maximum values of e c ,min = 0.0006, e C}max = 0.14, i.e., 
it changes from being almost perfectly circular to having 
an eccentricity of 0.14. The standard deviation figures 
are small, confirming that all 20 simulations display such 
oscillations. 

We have conducted a preliminary investigation of the 
TTVs expected in such a system, assuming that the sys- 
tem is edge-on to the line-of-sight. We find that over the 
course of 1-year of observations, TTVs of up to ±10 hours 
can occur. Such large variation would clearly be de- 



tectable, but would also have to be accurately accounted 
for in any initial search to ascertain whether the system 
transits at all (as any transit could easily be missed by 
many hours). 

We conclude that there are some dynamically fasci- 
nating Newtonian (interacting) interpretations of the ob- 
served radial velocity measurements, but that most so- 
lutions are unstable. This both adds to our suspicions 
that the two-companion model for the observed signal is 
incorrect and amplifies the extraordinary nature of the 
discovery if it is real. 

6. POTENTIAL FALSE ALARMS 

At this point in our analysis, we have an intriguing 
and exciting result: a clear, strong signal suggesting a 
super-planetary system in a near-perfect resonance with 
a strong likelihood of large and detectable dynamical ef- 
fects, and the potential for transits. 

But an application of a healthy skepticism reveals a few 
anomalies and coincidences tha t require considera tion. 

First, 3:1 resonances are rare (j Wright et alJiMl . Sec- 
ond, such resonances are only perfect on average on long 
time scales; the actual orbital elements tend to librate 
about the equilibrium solution, so at any given epoch 
the period commensuribility appears imperfect. These 
considerations make the apparent exactness of the com- 
mensuribility suspicious. Third, the fact that the data 
admit but a few dynamically stable solutions (only 1% of 
our simulations are stable, and even those are highly dy- 
namically active) suggests that we may not have settled 
on the correct solution. 

Further, the relative phases of the planets are suspi- 
ciously coincidental, with both planets exhibiting radial 
velocity extrema and zero-crossings simultaneously. To 
understand why this must be a coincidence, consider that 
for circular orbits these zero-crossings correspond to infe- 
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Fig. 7. — Orbital evolution for a sample long-term stable system. We select one of the long-term stable systems from our dynamical 
analysis and plot in detail the evolution of the orbital elements over 10 s years, plotting on the left-hand side the detailed evolution over a 
single year, in the center the stretch of data from 1-10 years, on the right the long term evolution from 10 years on to 10 5 years. In the 
bottom panel we plot the three resonant arguments 8i (red), 82 (gray) & 83 (blue, see Eqn[TJfor definitions), showing that at least one of 
the argument (82, gray) is librating, albeit with a relatively large amplitude. In the second row we plot the range of orbital distances q 
(pericenter, bottom of bars) a (semimajor axis, central line) and Q (apocenter, top of bars) for the b and c* components (black and red, 
respectively). The eccentricity of the more massive outer object changes very little, while the eccentricity of the c* component oscillates 
from ~ 0.01 — 0.1 over the course of ~ 0.5 years. Finally, the characteristic semi-major axis oscillations in a resonant system are present 
in this simulation (top row) and give rise to variations in the period ratio ~ 3% on ~ 0.5 year timescales (top plot). If the c* component 
is real this result raises the exciting possibility of detecting such period and eccentricity oscillations directly from RV measurements. 



rior and superior conjunction (i.e. transit and secondary 
eclipse for edge-on systems). As such, the zero-crossings 
of the b and c* components should only appear simulta- 
neous for a few special viewing angles. 

Finally, it seems inherently unlikely that an unusual 
system such as this, for which there is no good analog 
among the thousands of stars previously surveyed for 
planets, would be the first system to emerge from a new 
planet survey (although, as we will see in Section 16.81 
the selection criteria of MARVELS may make this not so 
unlikely as it seems). 

All of these oddities could be explained if the RV signal 
were not actually the sum of two Keplerians, but some 
other shape that was approximately fit by such a sum. 
Could the faint, AO companions to MARVELS-1 have 
something to do with this unusual result? 

Prudence requires both ruling out other explanations 
for the observed signal and understanding why MAR- 
VELS might be especially sensitive to such systems. In 
the face of the novelty of this system and the suspicions 
raised by the coincidences above, healthy skepticism de- 



mands that we rigorously consider all potential false pos- 
itives to confirm that each has truly been ruled out. 

6.1. Formal False Alarm Probability 

The formal false alarm probability for the planetary 
companion is virtually zero, due to its high amplitude 
compared to the r.m.s. residuals to the fit. We have 
confirmed this with a second -companion FAP c alcula- 
tion using bootstrapping (see iWright et al.ll2009l for a 
thorough description.) We subtracted the best- fit 6 d pe- 
riod Keplerian orbit from the velocities, and assumed the 
null hypothesis, that the remaining scatter in the resid- 
uals was unpatterned noise of some character (possibly 
non-Gaussian). We then redrew these residuals (with re- 
placement), added to them the original 6-d signal, and 
performed a thorough search of these artificial data for 
double-Keplerian solutions. We repeated this procedure 
to produce 1,000 trial data sets. We recorded the r.m.s. 
value for the residuals for each trial, and compared them 
with the residuals to the 2-Keplerian fit of the authentic 
data set. 

The sidelobes in the power spectrum of this data set 
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makes searching for the best 2-companion solution diffi- 
cult, since our L-M fitter (RVLIN) is easily trapped in 
local minima (i.e. sidelobes). We therefore thoroughly 
searched each bootstrapped dataset by testing the 10 
most significant peaks in the residual periodogram, suf- 
ficient to find the best fit in the unscrambled data. Even 
with this extra effort, trials in 1000 had an r.m.s. resid- 
ual as low as the actual set, consistent with the amplitude 
of the true signal. 

This result demonstrates that the 2-day signal is real 
in the sense that it is not a spurious periodogram peak 
introduced by the signal of the b component and unstruc- 
tured noise interacting with the window function of the 
observations. Demonstrating that the signal is due to a 
planet and not some spurious astrophysical, instrumen- 
tal, or methodological false positive requires eliminating 
alternative explanations (a task to which we devote the 
rest of this section) and verifying that the planet is phys- 
ically plausible (i.e. that the derived orbit is dynamically 
stable, which we did in § I5.3[) . 

6.2. Stellar Pulsations 

Practically speaking, any non-sinusoidal periodic ra- 
dial velocity signature measured at finite signal-to-noise 
ratio can be well fit by a finite sum of sinusoids. It is 
therefore prudent to consider whether the entire 2.6 km 
s _1 signal is in fact due to stellar oscillations, and the sig- 
nal from the apparent c* component is in fact merely the 
most important harmonic component of a fundamental 
mode, or even a pulsational mode of the star excited by 
resonant interaction with the brown dwarf companion. 

The lack of photometric variation (see Section|4]) would 
seem to rule out any pulsation mechanism as the ori- 
gin of the 2.6 km s _1 signal; nonetheless we have in- 
vestigated the expected photometric variation of many 
classes of variable stars, and conclude that the most im- 
portant false positive to consider for a star of this effec- 
tive temperature and periods near 1-10 d is that of a 
low-amplitude Cepheid (of either type), such as Pol aris. 

Polaris is well studied (see, e.g. JBruntt et al1l2008l for 
a summary). Its pulsational RV signature is sinusoidal, 
and its amplitude has decreased from 6 km s -1 to about 
0.5 km s _1 over the last century. Similarly, its photomet- 
ric amplitude has decreased from around 140 mmag to as 
low as 30 mmag in V band in the year 2000, an d as low 
as 10 mmag around 2005. Bruntt ct al. (20QM) find no 
evidence of ove rtones of the 4 d pul sation period in RVs 
or photometry. lUsenko et al.l (|2005l ) find that the pulsa- 
tional modes of Polaris are not obvious from the spectra 
in that the effective temperature of the star varies by no 
more than 100 K and does not corre l ate w ith pulsational 
phase (although lHatzes fc Cochran! {2000) find that the 
line bisector variations are detectable when the RV am- 
plitude is 1.5 km s _1 ). 

A low-amplitude Cepheid is thus a particularly insidi- 
ous false positive for a broad RV program, since it could 
in principle have a small (and easily overlooked) ^10 
mmag photometric variation, small line bisector varia- 
tion, essentially no T e g variation, and a large, easily de- 
tected ~ 1 km s _1 RV signal reminiscent of a hot Jupiter 
companion. 

There are several reasons this scenario, or a similar one, 
does not apply to MARVELS- 1. The first is that most 
secondary pulsation modes are not harmonics of the dom- 



inant pulsation mode, but occur at a period with some 
non-integer period ratio (beat Cepheids usually have pe- 
riod ratios near 0.7 or 0.8, not 1/3). The second is that 
our prototype false positive, Polaris, shows no overtones 
at any period ratio, integer o r otherwise . The third is 
that spectroscopic analysis of iLee et al.1 ()2011h conclu- 
sively shows that MARVELS- 1 is a dwarf, not a giant 
or supergiant, and so is not expected to have pulsational 
modes with periods near 6 d and in any case does not 
lie in the instability strip. Finally, our photometry in 
2] places upper limits on the variability of MARVELS- 1 
that are nearly two orders of magnitude more stringent 
than the minimum variability of Polaris, our prototype 
confounder. 

6.3. Stellar Activity and Star- Companion interactions 

While stellar activity can mimic the effects of planetary 
companions, it is an inadequate explanation in this case. 
MARVELS-1 does not show extraordinarily large Can H 
& K emission or large usini values, and so we do not 
expect variatio ns of order 150 m s -1 , as seen here for the 
c* component (| Wright! l2005h . and certainly not at the 
level of the signal of MARVELS-16. The stability of the 
photometry precludes large spot-induced errors, and the 
low v sin i indicates that the observed signals are likely 
at much shorter periods than the rotation period of the 
star. 

The 3:1 period commensuribility is also inconsistent 
with one of the signals being due to stochastic stellar 
activity, which naively should not occur at a harmonic 
of the orbital period of MARVELS- 16. 

That said, the large mass and close separation of 
MARVELS- 16 suggests that we consider interactions be- 
tween it and its host star as a possible source of unusual 
or novel effects on precise radial velocity and bisector 
measurements. Depending on the true mass and temper- 
ature of MARVELS-16, it may detectably heat the host 
star's surface, altering its spectrum and convective flows 
there; it could cause a prolate distortion of the host star; 
or it could perhaps resonantly excite a normally dormant 
pulsation mode in the star at 1.965 d. Any of these ef- 
fects could plausibly create a problematic RV signature 
with correlated bisector changes. 

However, once again the photometric stability of 
MARVELS-1, documented in 2] at all of these frequen- 
cies, makes these possibilities very unlikely. Any varia- 
tion in the output of MARVELS-1 of the sort described 
above strong enough to vary its spectral lines' positions 
by 100 ms _1 should also be apparent in the photometry, 
at least at the level of our photometric precision. 

Further, it is difficult to construct any such scenario 
that produces significant power at three times the fre- 
quency of the b component, but not significant power at 
twice the period (which, for instance, prolate distortions 
would generate). 

In short, any confounding effect along these lines would 
have to be previously undetected by radial velocity sur- 
veys, generate large (100 m s _1 ) RV variations at almost 
exactly three times the frequency of the b component 
(and at no other frequencies), and do so with no de- 
tectable photometric signature. We can think of no such 
effects. 

6.4. Is MARVELS-1 b a binary brown dwarf? 
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One source of false positive is that MARVEL S-l b 
is actually an equal mass binary brown dwarf with to- 
tal minimum mass 28 Mj U p, in a 3:2 orbit-orbit reso- 
nance. In this scenario, the variable quadropole moment 
of the binary brown dwarf would induce a small addi- 
tional Doppler signature on MARVELS- 1 at twice the 
frequency of the binary orbit, which we might observe as 
the 2-day signature. The binary would need to be nearly 
equal mass, or else the fundamental orbital frequency 
would also produce a detectable signature. 

Unfortunately, this particular false positive cannot be 
dispensed with quickly, since the parameters are rather 
well constrained and the predicted Doppler signature 
from such a scenario is similar to that seen in the ob- 
servations. We present here a quick demonstration of 
this fact, using a toy model. 

If we assume circular, unperturbed orbits for both the 
binary and the binary's orbit about MARVELS-1; denote 
the stellar mass M and the mass of each component of 
the binary m(= 14Mj U p); denote the semimajor axis 

of the orbit of the center of mass of the binary around 
MARVELS-1 a (corresponding to period P = 6 d); and 
denote the total orbital separation 2d of the binary com- 
ponents from each other; we can use Kepler's laws to 
derive a rough, order-of-magnitude expression for the ra- 
tio d/a: 



1 



9 M + 2m 



The semi-amplitude of the variation in the star's accel- 
eration due to this variable quadropole for small d/a can 
be approximated as 



~K 



rn 



9 M + 2m 



(2) 
(3) 



where K' represents the semiamplitude of the pertur- 
bative signal. For m/M ~ 4 x 1CV 2 (assuming a high 
inclination for which the total BD binary mass is actu- 
ally 80 Jupiter), we have K 1 / K ~ 0.03 or o ne part in 25. 
Our actual semiamplitude ratio from § 13.41 is one part in 
22, which is close enough that further analysis might be 
warranted. 

Our healthy skepticism thus finds significant purchase 
here. In order to claim that the c* component is a sep- 
arate planet and not a spurious residual from the signal 
of a binary brown dwarf one must perform a more de- 
tailed analysis and show that such a scenario is dynam- 
ically impossible. We did not perform such an analysis 
for this system here, because we have a more important 
confounding effect to consider. 

6.5. RV Signals From Fainter Stars 

An obvious potential source of problems is the presence 
of the two companion stars described in S I2.2.1[ whose 
spectra might contaminate that of MARVELS-1 and pro- 
duce complex and spurious velocity signatures. 

We might also suspect that the entire RV signature is 
due to pulsations of one of the background objects, of the 



sort discussed in H6.2I This explanation also suffers from 
the fact that the photometry in Sj4] detects no variability 
for this system; if the pulsation modes are so dilute as to 
be invisible in our photometry they should not appear as 
a 2.6 km s _1 signature in the RVs . This explanation also 
shares the difficulty that there are no classes of pulsing 
stars that exhibit the 3:1 period commensuribility we see. 

One general source of false positives with blended tar- 
gets is a faint stellar binary whose signal is diluted by 
a brighter coincident star. In this case the Doppler sig- 
nature of the binary is proportionally diluted, and the 
signal of the binary can be misinterpreted as a planetary 
signal (|Torres et al.ll2004bO . 

We can rule out the 0//9 companion, which lies near 
the edge of our 2" HRS fiber during nominal pointing. 
Seeing, tracking, and pointing variations should thus pro- 
duce variable levels of contamination, and thus a variable 
Doppler signature that would lack any particular period- 
icity. Also, the contamination should be at a completely 
different level in our Keck velocities (since HIRES is slit- 
fed). Since we see a consistent signal in both the HET 
and Keck velocities, if the signal is due to a binary, then 
it must be the binarity of the 0//15 object or another 
unresolved contaminant. 

The 3:1 period commensuribility further restricts the 
possibilities of a blend. If the source of the 100 m s _1 
amplitude signal is to be attributed to the binary mo- 
tion of one of the companion stars, we must accept the 
coincidence that the effect has a frequency surprisingly 
near a small integer multiple of the orbital frequency 
of MARVELS-16. Alternatively, if the entire RV sig- 
nature of MARVELS-1 is to be due to a contaminat- 
ing object, then we still must confront the novelty of 
this short-period 3:1 resonant system, except with sig- 
nificantly higher RV amplitudes (and thus masses and 
luminosities) for the objects (since their signal would be 
significantly diluted by MARVELS- 1 ) . 

In short, it is difficult to construct a plausible scenario 
whereby two signals in a near-perfect 3:1 resonance are 
the product of anything other than Doppler signals from 
the primary star, MARVELS-1. 

6.6. Effects of spectral contamination 

In this section, we consider other sources of spurious 
Doppler signatures from spectral contamination. In the 
discussion below we will imagine how a cross correla- 
tion procedure would be corrupted by the presence of a 
contaminating stellar spectrum. We will refer to a "tem- 
plate" spectrum (in practice, taken without the iodine 
cell in the beam) and an "epoch" spectrum (a subse- 
quent observation, through the iodine cell). The change 
in velocity measured for an "epoch" spectrum with re- 
spect to the "template" spectrum is measured as the ve- 
locity shift at which the cross-correlation function (CCF) 
is maximized. This peak occurs at Av — v e — vt , where 
v e is the velocity of the star at epoch, and v t was the 
velocity of the star during the template observation. 

If a companion star with radial velocity v c were con- 
taminating our spectra, then this contamination would 
produce a weak signal in the CCF at the velocity shift 
that aligns it with the template spectrum, Av = v c — Vt 
(this effect is weak both because the spectral types do 
not match and because the contaminating spectrum is 
weak.) There would also be a similar spurious peak 
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from the contamination in the template spectrum cor- 
relating with the primary spectrum in the epoch obser- 
vation, at Aw = v e — v c (there is also a second-order peak 
from the contaminating spectrum interacting with itself 
at Aw = 0.) 

In this scenario, as the spectrum of MARVELS-1 shifts 
due to the influence of the brown dwarf, its spectrum 
will periodically align with the contaminating spectrum. 
When this happens, the "true" peak in the CCF at 
At; = w e — v t will align with the spurious peak at v c — v t . 
At values of w e near w c , the two CCF peaks will be 
blended, and the peak of their summed shape, which 
will be asymmetric, will not be at v e — v t , but will rather 
be "pulled" towards v c — v t by an amount that depends 
on the amount of contamination and the magnitude of 

Vc - V t - 

Since the widths of these peaks are proportional to 
the line widths in the spectra, this effect will only be 
important for values of \v e — v c \ less than or similar to 
a typical line width, and the effect should decrease in 
magnitude rapidly for larger values. 

The second spurious peak (at Aw = v e — v c ) does not 
produce a similarly time-variable velocity anomaly; it 
will pull the true CCF peak (at Aw = v e — w<) by a con- 
stant amount, which will be subtracted as a zero-point 
offset in the resulting differential velocities. 

The time-variable component of this "peak pulling" 
could thus cause systematic errors in our measured ve- 
locities such that they are erroneously closer to the ve- 
locity of the contaminating star when the spectra are 
near alignment. This would manifest as an anomalous 
plateau near the preferred velocity, with a characteris- 
tic width determined by the typical line widths of the 
primary and contaminating lines. Since the amplitude of 
the RV variation from MARVELS-1 b is of similar magni- 
tude to the line widths of MARVELS-1, we would expect 
all of the velocities to be affected to some degree (except 
at the preferred velocity where the spectra are perfectly 
aligned), but most severe when the velocity difference is 
greatest. 

Although the AO companions to MARVELS-1 are too 
faint to contribute sufficient flux in the optical for this 
effect to be significant, we nonetheless explore it out of 
an abundance of caution. Below, we build a model for 
this effect and apply it to our Keck and HET velocities. 

6.6.1. A general model for "peak pulling" 

We have modeled such an anomaly as i> pcr t where the 
velocities show a systematic tendency towards a preferred 
velocity Vq when they are within a characteristic line 
width Aw of wo: 

w pert = fc(e-«^)/ A ^ 2 )(w -w) 

where k is the magnitude of the effect. 

We have generated artificial data assuming Aw = 3.5 
km s _1 , k = 0.25, and wo = as a perturbation on a 
sinusoidal signal of semiamplitude 3 km s , and per- 
formed a sinusoidal best fit to the resulting velocities; 
the result is presented in Figure HI Note the similarity of 
these curves to those in Figure |U 

Here, our healthy skepticism and abundant caution 
seems to have paid off. This "peak pulling" scenario nat- 
urally produces a residual signal to a sinusoidal fit whose 



TABLE 7 

Best Properties of the MARVELS-1 System 



Property Value Source 1 



Per (d) 5.895322 ± 0.0002 PPM 

T (JD-2440000) 15509.039 ± 0.07 PPM 

e 0.0002 to ou02 PPM 

w (°) 130 ± 0.4 PPM 

X Aa (km s- 1 ) 3.05t°'°i CCF 

K Ah (km s- 1 ) 3.mt° ° i CCF 

q l-20±goi CCF 

flux ratio ~ 0.2 CCF 

FWHM (Aa) (km s" 1 ) 7.26±g;g| CCF 

FWHM (Ab) (km s" 1 ) 7.31 ± 0.26 CCF 

i(°) 2.47 ± 0.04 §[6X6] 

A^Aa (M Q ) 1.25±0.06 §0 

M Ab (Mq) 1.04±0.05 M A Jq 

r.m.s. (m s 14.4 PPM 

Xl 0.97 PPM 



1 PPM: "Peak pulling model" from Section 16.6.11 
"CCF" : BIS and FWHM analysis of HET spectra from 
Section 16.6.41 

dominant mode is at exactly 3 times the frequency of the 
primary signal. 

This result can be understood from the symmetry pro- 
duced by a pulling near the velocity zero, which produces 
a signal that shares the symmetries of a sinusoid, but in 
this case is slightly taller and narrower. As Figure [8] illus- 
trates, the best-fit sinusoid and perturbed velocity curves 
meet at the w = points (because wo =0), and since the 
best fit will, by construction, follow the actual signal as 
closely as possible, it necessarily overestimates the true 
velocity near the w = points and underestimates it near 
the crests. The best-fit curve thus crosses the true curve 
twice between the nodes, producing six points with zero 
residuals per full period, and so an apparent 3:1 period 
commensuribility. 

The actual perturbative signal does not have this char- 
acteristic: the residuals to the true RV curve have a dif- 
ferent characteristic shape, as Figure [9] shows. 

We have fit the HET data (as reduced with the HET 
template) under this scenario by assuming that the data 
can be modeled as a Keplerian plus a w por t component 
with the free parameters given in Equation 16.6.11 The 
best-fit solution is superior (r.m.s. = 14 m s _1 ) to that 
given by a double-Keplerian to the HET data alone 
(r.m.s. = 19 m s _1 ), and is consistent with the errors. 
In this solution we find K = 3190 m s _1 , significantly 
higher than the fit value for K in a single-Keplerian 
solution, and an orbit consistent with circular (best fit 
e = 0.0002+ i ; I qq 2 ). The best fit perturbative parameters 
are wo = —20 m s _1 , k — 0.30, and Aw = 4.1 km s _1 , 
but are weakly constrained. 

6.6.2. What is the source of contamination? 

There are three obvious candidates for the origin of a 
contaminating spectrum: the 0'.'9 companion, the 0'.'15 
companion, and MARVELS-1 b. To first order, we ex- 
pect none of these to be responsible: the AO compan- 
ions are both 4 magnitudes fainter than MARVELS-1 in 
the NIR; if they are physically associated must be lower 
mass than MARVELS-1 (and so redder) and thus con- 
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Fig. 8. — Left: Velocities perturbed by a simple contamination model (black) and best circular fit (red). The blue dashed line illustrates 
the true, unperturbed velocity. Right: Residual velocities (perturbed minus fit) showing the nearly sinusoidal, triple-peaked residuals (i.e., 
black curve minus red curve from the left panel). The shape of these residuals are sensitive to the details of the contamination model, but 
will show a perfect 3:1 period commensuribility for special values, like the ones we have chosen here. We phase both diagrams to have 
phase zero at the velocity zero with positive slope of the primary component. Compare these figures with Figure [3] 




Phase Phase 

Fig. 9. — Left: The velocity perturbations v per t to the true RV curve from the same model used in Figure [8] (i.e. blue curve minus black 
curve from the left side of Figure l8l as opposed to the residuals to the fit RV curve, which are shown in the right panel of Figure [HJ. Right: 
Implied v pcI t curve from best fit to HET RVs of a Keplcrian-plus- "pulling"-model. This model fits the HET data significantly better than 
the double Keplerian model. We phase both diagrams to have phase zero at the velocity zero with positive slope of the primary component. 



tribute negligibly to the spectrum in the optical iodine 
region. MARVELS-1 b, being substellar, should have no 
appreciable optical emission; indeed, any object bound 
to MARVELS-1 must have a comparable luminosity to 
the star itself to be a spectral contaminant. Nonetheless, 
our healthy skepticism requires that we carefully exclude 
each possibility. 

The Keck velocities provide a powerful diagnostic, be- 
cause, unlike with the fiber-fed HRS, the apparent veloc- 
ity difference of MARVELS-1 and its contaminant will 
be a function of the position angle of slit during the 
measurement if the two objects have any measurable an- 
gular separation. Since the Keck 0'.'861 slit was always 
rotated to the pa rallactic angle during our observations 
( Filip penkol 119821 ). this would make the signal from the 
0'.'15 companion a strong function of hour angle of obser- 
vation. Since MARVELS-1 was observed at hour angles 
ranging from +0.6 to —3.1, this effect would make the 
velocity perturbations appear essentially random, with 
no relationship to the phase of MARVELS-1 b. Further, 
if the 0'.'9 companion were responsible, then its contribu- 
tion would be highly variable, since its light would often 
miss the narrow HIRES slit entirely. 



Interestingly, the Keck velocities fit the peak pulling 
model equally well (r.m.s. = 14 m s _1 ), although with a 
high reduced x 2 value due to the superior internal errors. 
The fact that the slit-fed HIRES sees a nearly identical 
signal to the fiber-fed HRS strongly suggests that the 
contaminating spectrum is located at small angular sep- 
aration, indeed within (f.'lb. 

Although the Keck velocities show a similar signal, fit- 
ting the Keck and HET data jointly does not produce a 
good fit in this model (r.m.s. = 33 m s _1 ), nor does fit- 
ting the HET data reduced with the Keck template alone 
(r.m.s. = 19 m s _1 ). This illustrates that the form of 
the velocity perturbations from the contaminating spec- 
trum is sensitive to the choice of template observation 
and the specific velocity reduction code used. This is 
perhaps not surprising, since the velocity perturbations 
are the result of a failure of the forward modeling process 
(which assumes uncontaminated spectra) and that dif- 
ferent templates will have the contaminating spectrum 
Doppler shifted with respect to the primary spectrum 
by different amounts (at the time of the Keck template, 
v = —2.0 km s _1 ; at the time of the HET template, 
v = +1.3 km s _1 ). 
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Finally, we consider the possibility that the spectrum 
of MARVELS- 1 b itself as a source of contamination, that 
is, that MARVELS-1 is a double-lined spectroscopic bi- 
nary, and our spurious signal at 2 d is a result of incorrect 
modeling of the spectrum as being only single-lined. In 
this case, the RV amplitude of the secondary spectrum 
could be quite large, depending on the mass ratio of the 
system. If this spectrum is contributing to Doppler sys- 
tematics, it might produce spurious (and small) peaks in 
the cross correlation function at its velocity during the 
template observation, at zero velocity phase (i.e. dur- 
ing superior and inferior conjunction, when its spectrum 
aligns with that of MARVELS-1), and perhaps also at 
the velocity of any other contaminating spectrum (e.g. 
from one of the AO companions) . 

The fact that we achieve different qualities of fit using 
the Keck template and the HET template is consistent 
with this picture, since the velocity of the companion 
spectrum in the template observation should be differ- 
ent between the two telescopes due to its large intrinsic 
motion. 

To check this possibility more rigorously, we con- 
structed a peak-pulling model where the peak is pulled 
toward a variable velocity consistent with that expected 
from MARVELS-1 b, instead of a constant velocity Vq. 
We find that the x 2 surface to fits to the HET data 
has a broad valley near its minimum, admitting signifi- 
cant motion of the secondary spectrum, proportional to 
Av. Indeed, we find that we cannot constrain the im- 
plied mass ratio at all, with equally good fits existing 
at (Av,m b /M) = (8.5km s"\ 0.65) and = (3.5 km s -1 , 
20.8). 

The former fit corresponds to a solution where 
MARVELS-1 b is stellar (that is, the system is a face- 
on binary, and the RV companion is not a brown-dwarf 
desert object, at all). The latter fit seems rather unphys- 
ical, but corresponds to an essentially fixed contaminant 
(i.e. it is equivalent to our original peak-pulling motion 
with no motion of the contaminating spectrum). This 
corresponds to a scenario in which there is a fifth object 
in the MARVELS-1, system, a chance alignment with a 
star of similar spectral type to MARVELS-1 but only 
15-30% of its flux, undetected in AO. 

To confirm one of these or similar unlikely scenarios, we 
must identify intruding lines themselves in the spectrum, 
and constrain whether they show significant Doppler mo- 
tion. To break the degeneracy we find in the peak pulling 
models, this analysis must involve an actual measure- 
ment of the line widths and bisectors in the observed 
spectra. 

6.6.3. Contamination Estimates and a Hunt for 
Contaminating Lines 

We first performed a perfunctory examination of the 
14 Keck spectra, which are of high signal-to-noise ratio 
and resolution, to find any obvious signs of peak pulling. 
Such an analysis would not necessarily be dispositive, be- 
cause the slit-fed HIRES line bisectors will be a strong 
function of the slit illumination profile, which may over- 
whelm any small line profile variations intrinsic to the 
source. Nonetheless, if the variations are sufficient to in- 
duce amplitudes of several hundred m s _1 , they may be 
visible. 

We examined two regions in the red (outside of the io- 



dine region): one region with strong telluric absorption 
and an adjacent order with several strong, deep stellar 
lines. We performed a simple cross correlation in pixel 
space in the telluric region to account for night-to-night 
lateral shifts in the spectral format, and interpolated the 
spectra onto a common scale. We then adopted an ap- 
proximate wavelength solution for the CCD (obtained 
from a ThAr exposure on another night). Next, we ap- 
plied the barycentric velocity shift to the wavelengths 
associated with each spectrum 

We selected three moderately deep and isolated stellar 
lines and fit Gaussians to their profiles, then inverted 
and normalized these profiles and interpolated them into 
velocity (log wavelength) space with zero velocity chosen 
to be the line center (as determined from the Gaussian 
fit). This approach allowed us to stack the line profiles 
for examination with improved signal-to-noise ratio. 

We show the stacked profiles for all 14 observations in 
Figure [TU1 with each component color-coded according to 
the measured radial velocity of MARVELS-1. The line 
profile is clearly undergoing variations, most especially 
in the line wings at the velocity extremes. This result 
is consistent with an underlying contaminating line at 
each position with velocity at nearly exactly the average 
of the line profile positions and with a slightly larger 
line width than that of MARVELS-1. In the frame of 
the figure (the frame of MARVELS-1) a contaminating 
spectrum protrudes in the wings at the extreme velocities 
(emerging most clearly in the wings near ±10 km s _1 ). 

6.6.4. CCF Bisector analysis 

Emboldened by the success of the peak pulling model 
at reproducing the observed RVs and its apparently suc- 
cessful prediction that there would be large line bisector 
variations, we proceed to perform a proper bisector anal- 
ysis of the fiber-fed HRS data, which does not suffer from 
slit illumination issues. The aim of the bisector analysis 
was to measure the large line profile variations predicted 
by the peak pulling model, and thereafter recreate the 
correlation between orbital phase and bisector change in 
an attempt to break the degeneracy between constant 
and variable velocity contaminants. In the context of 
planet-search surveys, there have been only a few in- 
stances where simulations of the bisector variation have 
been used t o convincingly dise ntangle blended radial ve- 
locities (e.g lSantos et al.ll2002l ). Conversely, the lack of a 
correlation between RV and bisector variations has some- 
times been used to constrain the parameters of a sys- 
tem or to rule out blend scenarios (jTorres et al.l l2004at 
IDi'az et al.ll20Tl . 

Line profile bisectors have conventionally been used to 
distinguish between true radial velocity signatures and 
spectral line asymmetries that m imic the existence of a 
companion (jToner fc Gravl [l988i). However, long expo- 
sure times or the stacking of many individual exposures 
are required in high-resolution spectroscopy to achieve 
the signal-to-noise ratio (SNR) required for reliable line 
bisector analysis, which makes it difficult to measure 
short-term variations. Instead, we scrutinize the bisec- 
tors of cross-correlation functions (CCFs) constructed 
from the cross-correlat ion of observed spec tra with nu- 
merical line masks fe.g. lBastiirk et aLll2011[ ). In addition 
to preserving line profile information, CCFs represent an 
"average" spectral line since they are based on multi- 
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Fig. 10. — Left: Stacked and normalized line profiles for 14 observations of MARVELS-1 from Keck. Colors indicate the measured radial 
velocity of MARVELS-1 (red indicating maximum redshift, etc.) Right: Detail in the line wings, showing the effect of the underlying 
contaminating spectrum. When the star is most redshifted (in the barycentric frame of the MARVELS-1 system, indicated by the red 
lines) there is a clear excess in the blue (negative velocity) wing of the line, and conversely on the red wing of the most blueshifted lines. 



pie lines of different elements (|Dall et al.l 120061) . Only- 
changes to the ensemble of spectral line profiles is re- 
flected in the shape of the CCF, making it less vulnerable 
to small-scale inconsistencies. 

In an effort to study the true CCF profile, we used ob- 
servations of MARVELS-1 on the HET HRS "red" CCD 
detector (6100 - 7600 A) which is beyond the iodine re- 
gion. The spectra consist of 25 echelle orders that were 
wavelength calibrated using the closest available Th-Ar 
hollow-cathode lamp calibration exposure taken during 
the night. Since the calibration exposures were neither si- 
multaneous nor bracketed, the wavelengths calculated by 
this method may suffer from instrumental shifts. How- 
ever, since we are only using these spectra for CCF bisec- 
tor analysis and not for precise radial velocity extraction, 
this lack of high precision in the wavelength calibration 
will not effect our analysis. Only the five most well- 
behaved echelle orders we used, in an attempt to mini- 
mize arbitrary bisector variability. We also rejected five 
epochs due to low signal to noise ratios (SNR), which 
produced faulty bisectors dominated by small-scale vari- 
ation. 

We computed the CCFs for each epoch by the cross- 
correlation of fully reduced and calibr ated spectra with 
a weighted 6?2 stella r template mask (|Pepe et al.l l2002t 
iBaranne et al.lfl996l) . The mask was created by us ex- 
pressly for this pur pose from an NSO FTS solar atlas 
(jKurucz et alJ[l984T ). Consisting of 230 lines spanning 
the wavelength region between 6000 - 7800 A, the mask 
has non-zero values coincident with the wavelengths of 
distinct unblended stellar lines, and reflects the relative 
depth of absorption lines against the local continuum. 
For the cross-correlation, we use all lines that have a 
depth of 5% or more with respect to the stellar contin- 
uum, since we lose lines by setting the mask to zero in 
regions where the stellar signature is overwhelmed by tel- 
luric lines. The width of the mask lines is adjustable, and 
set here to be 3 km s" 1 , based on the resolution setting 
of the spectrograph. 

Each spectral order was cross-correlated independently 
and the resulting CCFs added to attain a composite. The 
bisector was calculated as the loci of the midpoints of 
horizontal lines connecting the two wings of each com- 
posite CCF. Several measures have been used to quantify 



the shape of bisectors; we employ the bisector inverse 
slope (BIS), defined as the difference of average veloci- 
ties between 10-40% and 55-85% of the total CCF depth 
(jQueloz et al.l HoOlL The errors on the measurement of 
the BIS and the full-width at half-maximum (FWHM) 
were determined from the dispersion in each value be- 
tween different echelle orders. All of the measured BIS 
and FWHM values, along with their 1 — a uncertainties, 
are presented in Table [51 

The bulk motion of a star would cause bisectors to os- 
cillate around a mean bisector without changes in shape 
or orientation. However, the bisectors of MARVELS-1 
change dramatically in the period observed, producing 
a large range in the BIS fFig fTTj) . We also show ob- 
servations of a Dra in the same wavelength region and 
calibrated by the same method for comparison; these bi- 
sectors are extremely stable and show hardly any change. 
The contrast between these stars emphasizes the fact 
that information about the system is encoded into the 
pattern of the MARVELS-1 bisectors. 

Figure fL2T a) shows a conspicuous relation between the 
BIS (measured from the CCF) and the cosine of the or- 
bital phase angle, as calculated from the precise radial 
velocity results (extracted by the iodine method). Specif- 
ically, the measured velocity and orbital phase are related 
as v = k cos((f>), where k is a given velocity semiampli- 
tude. We choose to use the phase instead of directly 
using the RVs because the more precise iodine velocities 
are measured differently from the CCF method, where we 
simply measure the centroid of a fitted Gaussian. Even 
though the RVs produced by both methods are consis- 
tent within their respective uncertainties, we avoid any 
discrepancies by using a quantity that is agreed upon by 
both techniques. 

The small values of BIS near cos(0) = and their 
symmetry about that value indicates that the primary 
and contaminating spectra are aligned at v = 0. Figure 
[T2Tb) shows the FWHM values measured simultaneously 
with the BIS. This also has a minimum near cos(0) = 0, 
reaffirming that the spectra are aligned at v = 0. This 
is consistent with the result from our peak pulling model 
that vo ~ 0. 

6.6.5. A Two-Component CCF model 
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We attempted to reproduce this relation between BIS 
and cos(0) with a simple model. In keeping with the idea 
that spectral contamination from an unseen companion 
is leading to "peak pulling" , we simulated a scenario 
where a large Gaussian (Gi , representing the true CCF of 
MARVELS- 1 alone) shifts relative to a secondary Gaus- 
sian (G2, representing the hypothesized contaminant), 
and periodically aligns with it at some velocity shift. At 
each phase step we calculated the BIS and FWHM of 
the combined CCF, produced by the sum of G\ and G2. 
A similar simula t ion w ith two Gaussians was performed 
by iSantos et al.1 12002) to recreate a linear BIS-RV re- 



lationship. They were able to detect a potential brown 
dwarf around a member of a binary system, using a con- 
stant velocity G\ and four free parameters (including a 
variable G2). 

There are six free parameters in our model: the ratio 
of amplitudes of the two Gaussians (A 2 /Ai), the FWHM 
of Gi (wi), the FWHM of G 2 (w 2 ), the mass ratio (q = 
M2/M1), the velocity semiamplitude of G\ (Kx), and 
an adjustable offset between measured model BIS and 
observed BIS (Abis). This offset is included to allow for 
the fact that the observed CCF is not a perfect Gaussian, 
as assumed in the model. We explored the possibility 
of a constant velocity contaminant with this model, and 
found that it was not possible to successfully match both 
BIS and FWHM with this configuration, although the 
BIS could be reproduced by itself. Thus, the behavior of 
the line profiles rules out the presence of a fifth object in 
chance alignment with this system. 

A variable velocity contaminant renders an extremely 
viable model. In this scenario, the Gaussians undergo 
motion as in a binary system with a given mass ratio. 
In other words, the velocity semiamplitude of G2 is de- 
termined by the mass ratio and the motion of Gi, or 
K2 = —qK\. Figure ll"27 a) shows the result of our sim- 
ulation, where both BIS and FWHM are matched si- 
multaneously by mimicking a binary system. We arrive 
at our best fit param eters using an AMOEBA routine 
(jNelder fc Meadl 119651 ). and conclude that the observa- 
tions are modeled well by a binary system with q — 
1.20±g;g?, Kx = 3.05™ kms^ 1 , A 2 /A x = 18.93±°;1|%, 
wx = 7.26™ km s-\ w 2 = 7.31±g;|| km s" 1 , A B is = 
—0.11 ± 0.02 km s _1 . Errors quoted here are formal er- 
rors on the parameters, which do not include systematic 
effects or covariances with the other parameters. Note 
that A2/A1 here is a proxy for the flux ratio of the bi- 
nary. These formal uncertainties here were calculated by 
fixing all parameters, but one which was allowed to float, 
and determining the value of the floating parameter that 
increased x 2 by 1 above its minimum value (30.5, with 
34 degrees of freedom). 

The fact that we are able to model the line profile 
variability so convincingly suggests that the culprit for 
the observed anomaly in the MARVEL S-l RV residuals 
is in fact a massive, relatively bright, bound companion. 

6.6.6. Final best parameters for the MARVELS-1 binary 

system 

The mass ratio and amplitude of the RV primary com- 
ponent should be well represented by the q and Kx pa- 
rameters in our two-component CCF model. Interpreta- 
tion of A2/A1 as the flux ratio of the two stars is com- 



plicated by that quantity's dependence on the choice of 
lines used in the two models; since we have not extracted 
the spectrum of the secondary here we should anticipate 
that this is only a rough estimate of the true flux ratio. 

The widths of the Gaussians (w\,W2) should corre- 
spond to the actual FWHM of the lines in the individual 
spectra (these values were computed after subtraction of 
the instrumental profile), and would indicate an equa- 
torial rotational velocity v sin i = FWHM/V3 - 4 km 
s _1 if we assume that a rotation broadening kernel of a 
uniform intensity disk with no differential rotation dom- 
inates the line width. We have validated this estimate 
by applying our CCF code to synthetic spectra and the 
Solar spectrum convolved with rotational kernels or var- 
ious equatorial velocities to calibrate our method. We 
find that for the FWHM values we measure are consis- 
tent with ii sin i values of 3.5-4.5 km s _1 , favoring the 
higher values. 

The parameter Abis encapsulates information about 
the inherent bisector span of the blended spectrum when 
the two spectra are perfectly aligned in velocity space. In 
any case, the formal uncertainty on these model parame- 
ters underestimates the true uncertainties in the physical 
parameters of the system because of the crudeness of our 
model of stellar spectral CCFs as pure Gaussians. 

We calculate the inclination of the binary orbit from 
the binary mass function we infer from the CCF model: 



sin 3 i 



Kfq(l + q) 2 P 
2iGM Aa 



(4) 



where Ma 3 is the mass of the primary star. Propagation 
of errors using the formal uncertainty on q and K and 
assuming M Aa = 1.25 ± O.O6M (Section l2TT]l yields i 
2:47 ± 0?04, where the uncertainty in i is dominated by 
the uncertainty in the mass of the primary. 
We report our best system parameters in Table El 

6.7. Summary of Potential False Alarm Consideration 

We have found that the combination of the 3:1 period 
commcnsuribility of the signals, their relative and abso- 
lute strengths, and the stability of the photometry, rule 
out all scenarios that we have considered but three: a 
planetary companion in a 3:1 mean motion resonance, a 
binary brown dwarf in a 3:2 orbit-orbit resonance, and a 
15- 30% spectral contaminant with a systemic velocity 
nearly identical to that of MARVELS-1 from an unre- 
solved companion. 

Our line profile analysis definitively identifies 
MARVELS- 16 as the source of the contaminating 
spectrum, allowing us to dispose of the other hypothe- 
ses. The high inclination of this nearly- face-on binary 
dilutes the true Doppler signature of the orbit down to 
values consistent with a brown-dwarf-desert candidate, 
and the effects of the contaminating spectrum create 
residuals to a sinusoidal fit that mimic a 3:1 resonance. 

6.8. Similarities with Transit Searches 

How can we understand why such a rare face-on binary 
system happens to be the first system detected by MAR- 
VELS? Fortunately, we can achieve this understanding 
without a thorough analysis of the MARVELS selection 
effects, which is far beyond the scope of this work. 
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TABLE 8 

Observed values of BIS and FWHM for HET HRS data of MARVELS- 1 



time BIS UncertaintyBis FWHM UncertaintypwHM 



-2440000 (UTC) 


m s 1 


m s 1 


m s 


m s 


15177.6109 


0.06 


0.1 


8.97 


0.1 


15178.5997 


0.14 


0.1 


8.96 


0.1 


15180.7994 


0.15 


0.1 


9.00 


0.1 


15181.7885 


-0.01 


0.1 


9.03 


0.1 


15182.7878 


-0.37 


0.2 


9.40 


0.2 


15183.5788 


-0.06 


0.2 


8.99 


0.2 


15184.5782 


0.12 


0.1 


9.09 


0.1 


15185.5711 


0.75 


0.3 


9.23 


0.3 


15483.7415 


-0.20 


0.1 


9.21 


0.1 


15483.7528 


-0.22 


0.2 


9.25 


0.2 


15484.9464 


0.08 


0.1 


8.48 


0.1 


15484.9577 


0.07 


0.1 


8.68 


0.1 


15485.7397 


0.53 


0.2 


9.23 


0.2 


15485.7510 


0.56 


0.2 


9.26 


0.2 


15497.7059 


0.69 


0.3 


9.44 


0.3 


15498.7192 


0.57 


0.1 


9.27 


0.1 


15498.9215 


0.30 


0.2 


9.30 


0.2 


15500.6947 


-0.28 


0.1 


9.21 


0.1 


15500.9073 


-0.32 


0.2 


9.22 


0.2 


15501.6961 


-0.18 


0.1 


9.04 


0.1 


15507.6863 


-0.02 


0.2 


8.99 


0.2 


15510.6632 


0.30 


0.1 


9.20 


0.1 


15510.6715 


0.37 


0.1 


9.18 


0.1 


15522.6391 


0.15 


0.1 


9.06 


0.1 


15522.6469 


0.07 


0.1 


9.03 


0.1 


15522.6547 


0.24 


0.0 


8.85 


0.0 


15522.8492 


0.04 


0.1 


8.58 


0.1 


15522.8570 


0.02 


0.1 


8.70 


0.1 


15523.6459 


0.06 


0.1 


8.86 


0.1 


15524.6390 


-0.39 


0.2 


9.15 


0.2 


15527.6284 


0.57 


0.2 


9.14 


0.2 


15576.7008 


0.01 


0.2 


8.99 


0.2 


15577.7000 


-0.32 


0.1 


9.38 


0.1 


15777.9371 


-0.29 


0.1 


9.20 


0.1 


15779.9355 


0.08 


0.1 


8.92 


0.1 


15782.9385 


0.10 


0.1 


8.68 


0.1 


15783.9280 


-0.42 


0.2 


9.12 


0.2 


15790.9292 


0.10 


0.2 


8.71 


0.2 


15791.9106 


0.23 


0.1 


9.24 


0.1 


15792.9209 


0.70 


0.3 


9.28 


0.3 


15796.8955 


0.12 


0.1 


8.55 


0.1 



MARVELS has observed a total of 3,300 FGK stars 
with V=7.6-12 (a magnitude limited selection) in 2008- 
2012 with each observed about 27 times over 2 year win- 
dow. In contrast, other radial velocity planet searches 
have typically searched of order 1,000 stars. Because 
of its design, MARVELS necessarily targets stars whose 
multiplicity, variability, and spectra properties are less 
well-known beforehand. In addition, because the target 
stars are fainter, they are mostly likely to be blended 
with a significant contaminant (both because they are 
fainter and because they are more distant, and so there 
is a higher probability of a blended companion). MAR- 
VELS is thus especially sensitive to blend scenarios and 
binaries, relative to other RV planet searches, and needs 
to confirm its lowest amplitude detections in a manner 
similar to transit searches. 

In a similar vein, the MARVELS program has also re- 
cently discovered the double-lined highly eccentric spec- 
troscopic binary TYC-3010-1494-1, which masquerades 
as a modest-eccentricity brown dwarf candidate (Mack 
et al. 2013, AJ, accepted). 

Nonetheless, the probability of the particular blend 
scenario of MARVELS-1 seems to low for such an ob- 



ject to be expected in the MARVELS survey. A rough 
order of magnitude calculation would estimate a ~ 50% 
binary fraction for stars, that ~ 10% of binaries are close 
binaries, that perhaps ~ 10% have mass ratios such that 
the secondary would contaminate the spectrum, and that 
0.01% of systems would be within 3 degrees of face-on 
(so the RV amplitude of the orbit is comparable to the 
line widths). This would suggest that MARVELS should 
encounter 3300 x 0.5 x 0.1 x 0.1 x 0.0001 ~ 0.02 such sys- 
tems. Either these estimates are collectively off by two 
orders of magnitude, or the MARVELS survey was quite 
unlucky to have encountered this system. 

This situation is in some ways similar to that of transit 
searches, which have required extensive efforts to confirm 
and validate planet candidates, which must be sifted out 
from a large number of classes of false positive, in par- 
ticular blends with eclip sing binary stars. For instance, 
IMandushev et ahl {2005) found a similarly insidious form 
of false positive in the form of a eclipsing binary orbiting 
a slightly evolved F star that produced a nearly achro- 
matic photometric signal indistinguishable from a transit 
in their photometric wide-field transit survey, and pro- 
duced a similar "peak pulling" signal to the one we see 
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Fig. 11. — CCF bisectors, shifted by their CCF-measured radial velocities to align them for clarity. The colors are based on 
BIS values. Top: MARVELS-1. These bisectors change dramatically in both shape and orientation in the period observed. 
Bottom: a Dra. Notice the factor of ten difference in the velocity axis range, demonstrating the stability of the a Dra bisectors. 



here. 

7. NOMENCLATURE 

Nomenclature standards for multiple systems are sub- 
jective and tricky to apply consistently. If we attempt 
to follow the Washington Multiplicity Cata l og (W MC) 
standard as described in iRaghavan et al.l (|2010f ) and 
recommended by the International Astronomical Union 
(IAU), then we must still decide how to organize these 
(possibly unbound) companions hierarchically and in 
what order to assign component letters. 

The natural hierarchy here is unclear, because we 



might designate the 0'.'9 companion "B", and to group 
the other three components into a nearly-unresolved "A" 
system, or we might put the 0'.'15 and 0'.'9 compan- 
ions on equal footing and dub them "B" and "C" re- 
spectively. We choose the latter (consistent with our 
Figure [1]). This makes the face-on short period binary 
MARVELS-1 Aa and MARVELS-1 Ab (where the low- 
ercase letters follow the WMC convention for stars, but 
are conveniently consistent with the prior names for these 
objects, MARVELS-1 and MARVELS-1 6, respectively). 
We note that MARVELS-1 C is unlikely to be bound. 
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Fig. 12. — (a, left) Points showing the relation between ob- 
served CCF bisector (BIS) values and the cosine of the mea- 
sured phase angle (coscp) for MARVELS- 1 (for convenience, 
the phase here is defined such that phase zero (ip = 0) occurs 
at the maximum radial velocity for the best- fit circular orbit). 
The solid line shows the BIS calculated from the simple model 
consisting of two shifting Gaussians, using parameters cho- 
sen to fit the observed BIS and full-width at half-maximum 
(FWHM) values as a function of this orbital phase, (b, right) 
Points as in (a); line shows the full- width at half-maximum 
(FWHM) calculated simultaneously from the two-Gaussian 
model. In both plots the bisector inverse spans and widths of 
the model and observed CCFs were measured in an identical 
fashion. The error bars are determined from the measurement 
variation between echelle orders. 



8. CONCLUSIONS 

We have performed a thorough analysis of the complex 
MARVELS-1 system. Using adaptive optics imaging, we 
have identified a companion at 0'.'9 separation that ap- 
pears to be a foreground late M dwarf, and another at 
C//15 separation that appears to be a late K or early M 
dwarf associated with the primary, which we classify as 
a F9 dwarf. 

The primary was previously thought to host a short- 
period companion occupying the "brown-dwarf desert" , 
MARVELS-1 b, the first sub-stellar companion discov- 
ered with the MARVELS instrument. Follow-up ra- 
dial velocity measurements revealed strong deviations 
from a Keplerian solution, with amplitude ~ 100 m 
s _1 and a an orbital frequency exactly three times that 
of MARVELS-1 b. 



We have identified three extraordinary explanations for 
the observed radial velocity signature of MARVELS-1: 
a pair of substellar objects in a near-perfect 3:1 mean- 
motion resonance with strong dynamical interactions, a 
binary brown dwarf in a 3:2 resonance, and a blended 
stellar component contaminating the spectra, creating 
apparent residuals to a Keplerian solution at three times 
the observed period. 

Identification of strong line bisector variations consis- 
tent with a contaminating spectrum confirms that the 
final scenario as the correct one, and our detection of 
significant radial velocity motion of the contaminating 
spectrum confirms that it is due, in fact, to MARVELS-1 
b itself, which is actually a stellar binary companion in a 
face-on orbit. 

In this case, the concern raised over the unusual prop- 
erties of the system caught an unlikely and insidious form 
of spectral contamination. A routine check for line pro- 
file variations would have caught the problem early, but 
there is little motivation for such a system at Keck, where 
large line profiles variations are expected, even on RV 
stable stars, due to changes in the spectrograph and slit 
illumination. Bisector analysis is further hampered on 
unstabilized spectrographs by the lack of a precise wave- 
length scale outside the iodine region, and by the pres- 
ence of iodine lines within the iodine region. Only in 
cases such as MARVELS-1 with large (3 km s^ 1 ) RV 
variations, or with fiber-fed spectrographs, can bisector 
changes be obviously attributed to astrophysical, as op- 
posed to instrumental, effects. Nonetheless, such cases 
may be common among hot Jupiters discovered with Ke- 
pler, and in such cases a line bisector analysis from Keck 
spectra will be an important validation step. 

In the absence of such a check, the only indication 
that something was amiss was a "feeling" that the sys- 
tem was too unusual (the resonance was too perfect, the 
two-planet fit never good enough, the phases of the two 
"planets" too well matched) and the difficulty in identi- 
fying a good dynamical solution. 

In summary, MARVELS-1 appears to be a stellar 
triple, with one presumably bound companion seen in 
AO and the other detected by radial velocities and line 
bisector variations. A foreground apparent companion is 
separated from this system by 0'.'9 on the sky, bringing 
the total number of detected companions to the primary 
object to four. 
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